PythonによるscRNA-seq解析 その1¶

使用データ¶

マウス15.5日胚の眼球から取り出した網膜の細胞。10X Chromium 3′ v2で解析。2つのReplicatesでそれぞれ約2,600細胞。

Mouse developing retina from Lo Giudice et al. https://doi.org/10.1242/dev.178103

Lo Giudice, Q., Leleu, M., La Manno, G. & Fabre, P. J. Single-cell transcriptional logic of cell-fate specification and axon guidance in early-born retinal neurons. Development 146, dev178103 (2019).

Pythonに関係ない部分の解析¶

データのダウンロード¶

好きな方法で公共データベースからダウンロードする。レプリケイトがふたつあって、それぞれ以下のSRR IDに対応。

fastq-dump SRR8181428 --split-files --gzip
fastq-dump SRR8181429 --split-files --gzip

Cell Rangerによる解析¶

Cell Rangerでマウスゲノム（mm10）へのマッピング、遺伝子ごとのカウントを実行。

レプリケイトは（論文でそう名付けられてたので）E2, F2という名前のバッチに設定。

cellranger count --id=RetinalBatchE2 \
   --fastqs=/path/to/SRR8181428\
   --sample=SRR8181428 \
   --transcriptome=/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A

cellranger count --id=RetinalBatchF2 \
   --fastqs=/path/to/SRR8181429\
   --sample=SRR8181429 \
   --transcriptome=/path/to/Cell_Ranger_References/refdata-gex-mm10-2020-A

講習では、このCell Ranger解析の結果ディレクトリ（の一部、講習で使う部分だけを抜き出したもの）を配布している。

こんな感じの構成になっていて、Cell Rangerのアウトプットが見えるはず。velocytoは後半で解説。

In [1]:
!tree ./data

./data
├── RetinalBatchE2
│   ├── outs
│   │   ├── filtered_feature_bc_matrix
│   │   │   ├── barcodes.tsv.gz
│   │   │   ├── features.tsv.gz
│   │   │   └── matrix.mtx.gz
│   │   └── web_summary.html
│   └── velocyto
│       └── RetinalBatchE2.loom
├── RetinalBatchF2
│   ├── outs
│   │   ├── filtered_feature_bc_matrix
│   │   │   ├── barcodes.tsv.gz
│   │   │   ├── features.tsv.gz
│   │   │   └── matrix.mtx.gz
│   │   └── web_summary.html
│   └── velocyto
│       └── RetinalBatchF2.loom
├── retinal.h5ad
└── retinal_velo.loom

8 directories, 12 files

ライブラリのimport¶

定番のライブラリのimportに加えて、scanpyをscの別名でimportする。

Figのちょっとした設定。あと一応、筆者の環境における各ライブラリのバージョンは以下の通り。

In [2]:
import warnings
warnings.filterwarnings('ignore')

import os
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt
import seaborn as sns

sc.logging.print_header()
sc.settings.set_figure_params(dpi=100, facecolor='white')

scanpy==1.9.1 anndata==0.8.0 umap==0.5.3 numpy==1.22.3 scipy==1.9.3 pandas==1.5.2 scikit-learn==1.1.3 statsmodels==0.13.5 python-igraph==0.10.2 pynndescent==0.5.8

データの読み込み¶

Scanpyの特徴は、データフレームを拡張した *anndata* と呼ばれる オブジェクトを使う点にある。anndataを使うと、ひとつのオブジェクトに遺伝子発現量のデータ、サンプルや細胞のアノテーション、遺伝子の情報などをまとめて格納できる。 anndataを使うことで、実験の情報が詰まったひとつのオブジェクトに対して処理を次々に実行し、さらに処理結果をそこに蓄積していくことができる。

Scanpyでは、10x Genomicsのデータは、結果のディレクトリを指定することでそのままロードが可能な関数が用意されている。

ディレクトリには、3つのファイルが書き出されている。

ひとつは、個別の細胞を識別するバーコードが記述されたbarcodes.tsvという名前のテキストファイル。各行ごとにひとつのバーコードが記述されている。

ふたつめは、計測された遺伝子が記述されたgenes.tsvという名前のタブ区切りテキストファイル。このファイルは左側の列に遺伝子のENSEMBL Gene ID、右側の列の遺伝子のシンボルが記述されている。

最後にmatrix.mtxというテキストファイルが、各細胞（バーコード）について各遺伝子の発現がいくつ観測されたのかカウント情報をまとめたファイル。このファイルはMatrix Market Exchange Formatsという形式で記述されており、疎行列（含まれる値の多くがゼロであるような行列）を比較的コンパクトに記述するための形式となっている。

Scanpyの10xデータ用読み込み関数を使うと、これらを同時に読み込んで、適切に構造化されたオブジェクトができあがる。

In [3]:
adata_E2 = sc.read_10x_mtx(path='./data/RetinalBatchE2/outs/filtered_feature_bc_matrix/', cache=True)
adata_F2 = sc.read_10x_mtx(path='./data/RetinalBatchF2/outs/filtered_feature_bc_matrix/', cache=True)
In [4]:
adata_E2
Out[4]:
AnnData object with n_obs × n_vars = 3166 × 32285
    var: 'gene_ids', 'feature_types'
In [5]:
adata_F2
Out[5]:
AnnData object with n_obs × n_vars = 3299 × 32285
    var: 'gene_ids', 'feature_types'

2つのバッチをひとつのanndataオブジェクトに統合する。バッチのラベルもここで設定。

In [6]:
adata = adata_E2.concatenate(adata_F2, batch_categories=['E2', 'F2'])

観測値（細胞）に関するデータは、obs でアクセスする。 observationsの略。このデータフレームに細胞に関するデータを追加していくことができる。

In [7]:
adata.obs
Out[7]:
batch
AAACCTGAGATGTCGG-1-E2 E2
AAACCTGCAATCCAAC-1-E2 E2
AAACCTGCACATGGGA-1-E2 E2
AAACCTGCAGCAGTTT-1-E2 E2
AAACCTGGTTCCTCCA-1-E2 E2
... ...
TTTGTCACATCGATGT-1-F2 F2
TTTGTCAGTAGAGGAA-1-F2 F2
TTTGTCAGTCATCCCT-1-F2 F2
TTTGTCATCAGTTCGA-1-F2 F2
TTTGTCATCCGAATGT-1-F2 F2

6465 rows × 1 columns

変数（遺伝子）に関するデータは、var でアクセスする。variable の略。現在は、遺伝子のシンボル（名前）と、Ensembl Gene IDなどが登録されている。

In [8]:
adata.var
Out[8]:
gene_ids feature_types
Xkr4 ENSMUSG00000051951 Gene Expression
Gm1992 ENSMUSG00000089699 Gene Expression
Gm19938 ENSMUSG00000102331 Gene Expression
Gm37381 ENSMUSG00000102343 Gene Expression
Rp1 ENSMUSG00000025900 Gene Expression
... ... ...
AC124606.1 ENSMUSG00000095523 Gene Expression
AC133095.2 ENSMUSG00000095475 Gene Expression
AC133095.1 ENSMUSG00000094855 Gene Expression
AC234645.1 ENSMUSG00000095019 Gene Expression
AC149090.1 ENSMUSG00000095041 Gene Expression

32285 rows × 2 columns

カウントテーブルには、以下のように *X* でアクセスできる。疎行列として格納されている。

In [9]:
adata.X
Out[9]:
<6465x32285 sparse matrix of type '<class 'numpy.float32'>'
	with 10273670 stored elements in Compressed Sparse Row format>
In [10]:
adata
Out[10]:
AnnData object with n_obs × n_vars = 6465 × 32285
    obs: 'batch'
    var: 'gene_ids', 'feature_types'

前処理¶

観測された遺伝子が極端に少ない細胞、割り当てられた細胞が極端に少ない遺伝子をフィルタリング。

In [11]:
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=20)
In [12]:
adata
Out[12]:
AnnData object with n_obs × n_vars = 6462 × 13332
    obs: 'batch', 'n_genes'
    var: 'gene_ids', 'feature_types', 'n_cells'

ミトコンドリア遺伝子の発現割合でフィルタリングする予定なので、遺伝子側のメタデータにミトコンドリア遺伝子か否かの情報を付与。

In [13]:
adata.var['mt'] = adata.var_names.str.startswith('mt-')

以下は、クオリティコントロールに関連したいくつかの指標を一挙に計算してくれる便利な関数。

qc_varsに遺伝子側のメタデータラベルを設定すると、それについてもカウントとパーセンテージを勝手に計算してくれる。

In [14]:
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
In [15]:
adata.obs
Out[15]:
batch n_genes n_genes_by_counts total_counts total_counts_mt pct_counts_mt
AAACCTGAGATGTCGG-1-E2 E2 1432 1427 3263.0 111.0 3.401778
AAACCTGCAATCCAAC-1-E2 E2 1297 1292 2568.0 97.0 3.777258
AAACCTGCACATGGGA-1-E2 E2 261 261 572.0 47.0 8.216784
AAACCTGCAGCAGTTT-1-E2 E2 360 359 624.0 21.0 3.365385
AAACCTGGTTCCTCCA-1-E2 E2 3072 3057 10505.0 283.0 2.693955
... ... ... ... ... ... ...
TTTGTCACATCGATGT-1-F2 F2 1916 1908 3714.0 136.0 3.661820
TTTGTCAGTAGAGGAA-1-F2 F2 2010 2007 5526.0 176.0 3.184944
TTTGTCAGTCATCCCT-1-F2 F2 1760 1757 3718.0 121.0 3.254438
TTTGTCATCAGTTCGA-1-F2 F2 1175 1170 2064.0 54.0 2.616279
TTTGTCATCCGAATGT-1-F2 F2 443 441 865.0 305.0 35.260117

6462 rows × 6 columns

それぞれの分布をプロット。

In [16]:
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True)
In [17]:
sc.pl.scatter(adata, 'total_counts', 'n_genes_by_counts', color='pct_counts_mt', size=40)
In [18]:
fig = plt.figure()
sns.displot(adata.obs['pct_counts_mt'][adata.obs['pct_counts_mt'] < 10], kde=False)
plt.show()

<Figure size 400x400 with 0 Axes>

分布を見て方針を決めて、フィルタリングを実行。

In [19]:
print('Total number of cells: {:d}'.format(adata.n_obs))

sc.pp.filter_cells(adata, min_counts = 2000)
print('Number of cells after min count filter: {:d}'.format(adata.n_obs))

sc.pp.filter_cells(adata, max_counts = 13000)
print('Number of cells after max count filter: {:d}'.format(adata.n_obs))

sc.pp.filter_cells(adata, min_genes = 1000)
print('Number of cells after gene filter: {:d}'.format(adata.n_obs))

adata = adata[adata.obs['pct_counts_mt'] < 6]
print('Number of cells after MT filter: {:d}'.format(adata.n_obs))

Total number of cells: 6462
Number of cells after min count filter: 4728
Number of cells after max count filter: 4706
Number of cells after gene filter: 4686
Number of cells after MT filter: 4563
In [20]:
adata
Out[20]:
View of AnnData object with n_obs × n_vars = 4563 × 13332
    obs: 'batch', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'n_counts'
    var: 'gene_ids', 'feature_types', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'

正規化や対数変換などの前に、この段階のカウントマトリックスを別のレイヤーに退避させておく。あとで確率モデル推定のときに必要になる。レイヤーに保管したデータは基本的に以降の操作の影響を受けないが、細胞や遺伝子のslicingは適用されるので注意。（常にマトリックスのshapeは一致する）

In [21]:
adata.layers['counts'] = adata.X.copy()

正規化¶

ライブラリサイズによる正規化、対数変換などの前処理は、scanpy.pp以下にいくつか便利な関数がある。

ここでは、細胞ごとのカウントの和が10,000になるように正規化してから、対数変換。

In [22]:
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)
sc.pp.log1p(adata)

この段階のデータ（いろいろとややこしい変換をしていない「ピュア」なデータ）も、あとでプロットや結果解釈のときに使いたいので退避させておく。rawは特別なレイヤーで、数値だけでなくobsやvarのメタデータも含めてフリーズしてくれる。また遺伝子側のsliceの影響を受けない（細胞側のsliceは適用される）

In [23]:
adata.raw = adata

特徴量選択（発現量の変動が大きい遺伝子）¶

発現量変動の大きい遺伝子のみを抽出して、データのサイズを小さくする。

In [24]:
sc.pp.highly_variable_genes(adata, flavor='seurat_v3', n_top_genes=2000)
print('\n','Number of highly variable genes: {:d}'.format(np.sum(adata.var['highly_variable'])))

 Number of highly variable genes: 2000
In [25]:
sc.pl.highly_variable_genes(adata)

計算結果は自動的に adata.var 、つまり、遺伝子に関するメタデータを格納したオブジェクトに追加される。 highly_variable の項目が True の遺伝子が高発現変動遺伝子。

In [26]:
adata.var
Out[26]:
gene_ids feature_types n_cells mt n_cells_by_counts mean_counts pct_dropout_by_counts total_counts highly_variable highly_variable_rank means variances variances_norm
Xkr4 ENSMUSG00000051951 Gene Expression 89 False 89 0.013773 98.622717 89.0 False NaN 0.013526 0.013121 0.834873
Gm19938 ENSMUSG00000102331 Gene Expression 578 False 578 0.103528 91.055401 669.0 False NaN 0.118130 0.129846 0.984716
Mrpl15 ENSMUSG00000033845 Gene Expression 1899 False 1899 0.375271 70.612813 2425.0 False NaN 0.471746 0.418216 0.943994
Lypla1 ENSMUSG00000025903 Gene Expression 1284 False 1284 0.241566 80.129991 1561.0 False NaN 0.302353 0.307116 0.989431
Tcea1 ENSMUSG00000033813 Gene Expression 2614 False 2614 0.631848 59.548128 4083.0 False NaN 0.704444 0.549619 0.948720
... ... ... ... ... ... ... ... ... ... ... ... ... ...
CAAA01118383.1 ENSMUSG00000063897 Gene Expression 96 False 96 0.015166 98.514392 98.0 False NaN 0.021325 0.025196 1.020294
Vamp7 ENSMUSG00000051412 Gene Expression 1606 False 1606 0.314144 75.147013 2030.0 False NaN 0.395332 0.365986 0.947294
Tmlhe ENSMUSG00000079834 Gene Expression 46 False 46 0.007119 99.288146 46.0 False NaN 0.007515 0.008047 0.921716
4933409K07Rik ENSMUSG00000095552 Gene Expression 20 False 20 0.003095 99.690498 20.0 False NaN 0.003088 0.002459 0.797902
AC149090.1 ENSMUSG00000095041 Gene Expression 1960 False 1960 0.475085 69.668833 3070.0 False NaN 0.462647 0.455644 1.044050

13332 rows × 13 columns

主成分分析（PCA）¶

PCAはscanpyの前処理関数で簡単に実行できる。とりあえず、50次元まで落としてみる。 use_highly_variableのフラグをオンにすると、遺伝子全体ではなく前項で決定した高発現変動遺伝子のみを使って次元削減をする。

In [27]:
sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack')

PCAで落とした座標は、観察値のメタデータを格納する *obsm* に自動的に入る。obsm は Multi-dimensional annotation of observations の略。複数の次元でひとつの何かを表現するような観測値のメタデータがここに格納される。

In [28]:
print(adata.obsm['X_pca'].shape)

(4563, 50)

プロットしてみる。Scanpyでは基本的に「前処理」（Preprocessing）に関わる関数がscanpy.ppに、「プロット」（Plot）に関わる関数がscanpy.plに入っている。

In [29]:
sc.pl.pca(adata, color='total_counts')
In [30]:
sc.pl.pca(adata, color='total_counts', 
          components=['1,2', '2,3', '1,3'])

t分布型確率的近傍埋め込み（t-SNE）¶

まず、t-SNE, UMAP共通のステップとして、データから「近傍グラフ」（neighborhood graph）の構築が必要。

In [31]:
sc.pp.neighbors(adata)

データ点間の接続関係（細胞の近傍関係）は、全細胞vs.全細胞のペアの情報を記録する *obsp* （Pairwise annotation of observations の略）に格納される。

In [32]:
adata.obsp['connectivities']
Out[32]:
<4563x4563 sparse matrix of type '<class 'numpy.float32'>'
	with 94424 stored elements in Compressed Sparse Row format>

接続関係を元にしてt-SNEを実行。

In [33]:
sc.tl.tsne(adata)

t-SNE結果のプロット。いまのところ特に色つけるようなデータも無いので、適当に細胞ごとのカウントで色分けをした。遺伝子名をcolorに指定するとその遺伝子の発現量によるプロットも可能。

In [34]:
sc.pl.tsne(adata, color='total_counts')

UMAP¶

近傍グラフはすでに計算しているので、scanpy.tlのumap関数を使えばオーケー。

In [35]:
sc.tl.umap(adata)

プロットはt-SNEとほとんど同じ。scanpy.pl以下にUMAP用の描画関数がある。

In [36]:
sc.pl.umap(adata, color='total_counts')

あきらかに似た形状の、ふたつの構造がある。。。

各細胞が由来するバッチで色分けしてみると・・・

In [37]:
sc.pl.umap(adata, color='batch')

ということで、このふたつの構造はもろにバッチ効果の結果だった。なんらかの方法で補正が必要だが、それは後述。

クラスタリング¶

Leidenクラスタリングを実行してみる。モジュラリティの計算に影響を与える "resolution" パラメータが存在する。

In [38]:
sc.tl.leiden(adata, resolution=0.5, key_added='leiden_r0.5')

クラスタリングの結果は観測値のメタデータ（obs）の中に格納される。

In [39]:
adata.obs
Out[39]:
batch n_genes n_genes_by_counts total_counts total_counts_mt pct_counts_mt n_counts leiden_r0.5
AAACCTGAGATGTCGG-1-E2 E2 1427 1427 3263.0 111.0 3.401778 3263.0 1
AAACCTGCAATCCAAC-1-E2 E2 1292 1292 2568.0 97.0 3.777258 2568.0 1
AAACCTGGTTCCTCCA-1-E2 E2 3057 3057 10505.0 283.0 2.693955 10505.0 2
AAACCTGTCCAATGGT-1-E2 E2 1449 1449 2781.0 91.0 3.272204 2781.0 7
AAACGGGAGGCAATTA-1-E2 E2 2772 2772 9240.0 224.0 2.424242 9240.0 1
... ... ... ... ... ... ... ... ...
TTTGTCAAGGGCTCTC-1-F2 F2 2393 2393 6209.0 143.0 2.303108 6209.0 5
TTTGTCACATCGATGT-1-F2 F2 1908 1908 3714.0 136.0 3.661820 3714.0 6
TTTGTCAGTAGAGGAA-1-F2 F2 2007 2007 5526.0 176.0 3.184944 5526.0 0
TTTGTCAGTCATCCCT-1-F2 F2 1757 1757 3718.0 121.0 3.254438 3718.0 5
TTTGTCATCAGTTCGA-1-F2 F2 1170 1170 2064.0 54.0 2.616279 2064.0 6

4563 rows × 8 columns

プロットする点の座標じたいはUMAPの結果なので、プロットはUMAP版の関数を使い、色分けだけをクラスタリング結果で指定すればいい。

In [40]:
sc.pl.umap(adata,
           color=['batch', 'leiden_r0.5'],
           ncols=2,
           frameon=False)

深層生成モデルの利用¶

scVIが仮定しているデータの生成プロセスは以下の確率モデルに基づく。

全細胞数がN, 細胞のインデックスが $n \in \{1..N\}$、全遺伝子数がG、遺伝子のインデックスが $g \in \{1..G\}$、細胞nのバッチIDが$s_n$であるとしたとき、細胞nの遺伝子gのカウントは、

$$z_n \sim Normal(0, I)$$$$l_n \sim logNormal(l_\mu, l_\sigma^2)$$$$\rho_n = f_w (z_n, s_n)$$$$w_{ng} = Gamma(\rho_n^g, \theta)$$$$y_{ng} = Poisson(l_n w_{ng})$$$$h_{ng} = Bernoulli( f_h^g(z_n, s_n) )$$$$ \begin{equation} x_{ng} = \begin{cases} y_{ng} & \text{if}\ h_{ng}=0 \\ 0 & \text{otherwise} \end{cases} \end{equation} $$

$l_\mu, l_\sigma$は、スケーリングファクターでバッチの数をBとしたとき$l_\mu,l_\sigma \in \mathbb{R}^B$で、実際のバッチごとの平均と分散に設定しておく。

ガンマ分布のパラメータ、ベルヌーイ分布のパラメータは、正規分布からサンプリングした$z_n$を使ってニューラルネットワークで表現する（VAEにおけるreparametrization trick）。

最終的なカウントはガンマとポアソンの複合分布なので、負の二項分布でモデル化してることになる。さらにベルヌーイ分布で、scRNA-seqにありがちなdrop-outイベントを表現。これらの組み合わせで、ゼロ過剰負の二項分布(Zero-inflated negative binomial distribution; ZINB)を表現している。

事後分布はVAEをSGDで最適化して変分推論。

学習後、$z_n$にアクセスすれば、バッチの効果が除去された細胞の潜在ベクトルが得られるし、$\rho_{ng}$にアクセスすればdrop-outなどのイベントが生じなかった場合の遺伝子発現の期待値が得られる。

In [41]:
import scvi

Global seed set to 0
/root/anaconda3/envs/pags/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:53: LightningDeprecationWarning: pytorch_lightning.utilities.warnings.rank_zero_deprecation has been deprecated in v1.6 and will be removed in v1.8. Use the equivalent function from the pytorch_lightning.utilities.rank_zero module instead.
  new_rank_zero_deprecation(
/root/anaconda3/envs/pags/lib/python3.9/site-packages/pytorch_lightning/utilities/warnings.py:58: LightningDeprecationWarning: The `pytorch_lightning.loggers.base.rank_zero_experiment` is deprecated in v1.7 and will be removed in v1.9. Please use `pytorch_lightning.loggers.logger.rank_zero_experiment` instead.
  return new_rank_zero_deprecation(*args, **kwargs)

データを格納したanndataオブジェクトと、バッチのIDを指定したobsのカラムのラベルを指定して、モデルをセットアップする。カウントデータの確率モデルなので、（正規化や対数変換したマトリックスではなく）カウントデータを格納したレイヤーを指定する必要があることに注意。

In [42]:
scvi.model.SCVI.setup_anndata(
    adata,
    layer='counts',
    batch_key='batch',
)

WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)

モデルを作る。ここで中間層のノード数とか潜在ベクトルの次元サイズなど設定できるが、とくに変更しなくていいと思う。

In [43]:
model = scvi.model.SCVI(adata)

scVIモデルのトレーニング（バッチ補正）¶

学習を実行する。GPU使わない場合はめっちゃ時間かかるので今回の講習では割愛。学習済みのモデルパラメータを用意したのでそれをロードして学習できたことにする。

In [44]:
# model.train()
In [45]:
# model.save('./models/scVI_model', overwrite=True)
In [46]:
model = scvi.model.SCVI.load('./models/scVI_model', adata=adata)

INFO     File ./models/scVI_model/model.pt already downloaded

学習された細胞ごとの潜在表現$z_n$は以下の関数で取得できる。PCAやUMAPの座標と同じように、obsmに格納しておく。

In [47]:
adata.obsm['X_scVI'] = model.get_latent_representation()

発現量期待値 $\rho_{ng}$ 、つまりdenoiseされた発現量テーブルは以下の関数で取得できる。適当なライブラリサイズで規格化してくれる。この数値はあとで使いたいので別のレイヤーに格納しておく。

In [48]:
adata.layers['scvi_normalized'] = model.get_normalized_expression(library_size=1e4)
In [49]:
adata
Out[49]:
AnnData object with n_obs × n_vars = 4563 × 13332
    obs: 'batch', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'n_counts', 'leiden_r0.5', '_scvi_batch', '_scvi_labels'
    var: 'gene_ids', 'feature_types', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm'
    uns: 'log1p', 'hvg', 'pca', 'neighbors', 'tsne', 'umap', 'batch_colors', 'leiden', 'leiden_r0.5_colors', '_scvi_uuid', '_scvi_manager_uuid'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_scVI'
    varm: 'PCs'
    layers: 'counts', 'scvi_normalized'
    obsp: 'distances', 'connectivities'

scVI潜在表現+UMAPによる次元削減¶

潜在表現でちゃんとバッチ効果が補正されたかどうか確認してみる。

近傍グラフ計算の scanpy.pp.neighbors は、デフォルトではPCAで計算した主成分（adata.obsm['X_pca']）を元にグラフを構築する。ここではPCA結果の座標ではなく、scVIで推定された細胞の潜在ベクトルを指定してグラフ構築してみる。バッチ効果が補正されたグラフが構築されるはず。

In [50]:
sc.pp.neighbors(adata,
                n_neighbors=30,
                use_rep="X_scVI")
sc.tl.umap(adata, min_dist=0.5)
sc.pl.umap(adata, color='batch')

scVI潜在表現+Leidenによるクラスタリング¶

同じように、scVI潜在表現で構成された近傍グラフを元にLeidenクラスタリングを実行してみる。上のセルで scanpy.pp.neighbors を実行したことにより、adata.obspにはscVI潜在表現で構築された近傍グラフが入っているので、ここでは普通に scanpy.tl.leiden を実行すればいい。

In [51]:
sc.tl.leiden(adata, key_added="leiden_scVI", resolution=1.0)
sc.pl.umap(adata,
           color=['batch', 'leiden_scVI'],
           ncols=2,
           frameon=False)

Doublet検出¶

SOLOによるDoublet検出。

Bernstein, Nicholas J., et al. "Solo: doublet identification in single-cell RNA-Seq via semi-supervised deep learning." Cell systems 11.1 (2020): 95-101.

時間に余裕があればやる。

In [52]:
results = []
for batch in ['F2', 'E2']:
    tmp_solo_model = scvi.external.SOLO.from_scvi_model(
        model, restrict_to_batch=batch)
    tmp_solo_model.train()
    result = tmp_solo_model.predict(soft=False)
    # 意図は不明だがSOLO予測のindexになぜか "-0" が付加されるので消しておく
    result.index = result.index.str.replace("-0$", "", regex=True)
    results.append(result)
results = pd.concat(results)

INFO     Creating doublets, preparing SOLO model.                                            

GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs

Epoch 79/400:  20%|█▉        | 79/400 [00:13<00:54,  5.90it/s, loss=0.144, v_num=1]
Monitored metric validation_loss did not improve in the last 30 records. Best score: 0.176. Signaling Trainer to stop.
INFO     Creating doublets, preparing SOLO model.                                            

GPU available: False, used: False
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs

Epoch 104/400:  26%|██▌       | 104/400 [00:15<00:45,  6.56it/s, loss=0.131, v_num=1]
Monitored metric validation_loss did not improve in the last 30 records. Best score: 0.158. Signaling Trainer to stop.
In [53]:
adata.obs['SOLO_prediction'] = results[adata.obs.index]
In [54]:
adata.obs['SOLO_prediction'].value_counts()
Out[54]:
singlet    4109
doublet     454
Name: SOLO_prediction, dtype: int64
In [55]:
sc.pl.umap(adata, color='SOLO_prediction', frameon=False)

DEG解析¶

仮説検定ではなく、ベイズ統計の枠組みでモデル比較を行う。基本的にはある集団と別の集団それぞれの遺伝子発現量期待値（VAEでdenoiseした発現量$\rho_{ng}$）のlog fold changeの値が、ある閾値以下となるか、それ以上の変化があるか、のふたつのモデルで比較する。

groupbyにラベルを指定すると、そのクラス（この場合はleidenクラスタリングのひとつのクラスタ）とその他全部で自動的に比較してくれる。比較のグループを一部に制限したり、ペアを指定したりいろいろと複雑に組み合わせは指定できる。詳細はAPIのドキュメントに。

In [56]:
DEs = model.differential_expression(groupby='leiden_scVI')
DEs.head()

DE...: 100%|██████████| 11/11 [00:54<00:00,  4.95s/it]
Out[56]:
proba_de proba_not_de bayes_factor scale1 scale2 pseudocounts delta lfc_mean lfc_median lfc_std ... raw_mean1 raw_mean2 non_zeros_proportion1 non_zeros_proportion2 raw_normalized_mean1 raw_normalized_mean2 is_de_fdr_0.05 comparison group1 group2
Carmil2 0.9836 0.0164 4.093937 2.658159e-07 0.000005 0.0 0.25 -5.505758 -5.591664 5.274880 ... 0.001264 0.015907 0.001264 0.015642 0.001908 0.021632 True 0 vs Rest 0 Rest
Pcdhb18 0.9822 0.0178 4.010596 8.968271e-08 0.000006 0.0 0.25 -6.042808 -6.083243 6.751517 ... 0.000000 0.006628 0.000000 0.006098 0.000000 0.012190 True 0 vs Rest 0 Rest
Cntn5 0.9812 0.0188 3.954919 2.587765e-07 0.000017 0.0 0.25 -5.148154 -5.396956 4.268237 ... 0.001264 0.074231 0.001264 0.056469 0.002287 0.091536 True 0 vs Rest 0 Rest
Gfra3 0.9808 0.0192 3.933458 2.856224e-08 0.000010 0.0 0.25 -6.871420 -7.264495 5.107664 ... 0.000000 0.022269 0.000000 0.016702 0.000000 0.037028 True 0 vs Rest 0 Rest
Psd2 0.9808 0.0192 3.933458 9.740267e-08 0.000005 0.0 0.25 -6.067220 -6.589231 5.682239 ... 0.000000 0.016437 0.000000 0.015642 0.000000 0.023997 True 0 vs Rest 0 Rest

5 rows × 22 columns

それぞれのクラスタ「特異的」遺伝子、DEG確率の高い順に上から3つくらい見てみる。

In [57]:
markers = {}
for i, c in enumerate(adata.obs['leiden_scVI'].unique()):
    comparison_label = "{} vs Rest".format(c)
    cluster_df = DEs.loc[DEs['comparison'] == comparison_label]
    cluster_df = cluster_df[cluster_df['lfc_mean'] > 0]
    cluster_df = cluster_df[cluster_df['bayes_factor'] > 3]
    cluster_df = cluster_df[cluster_df['non_zeros_proportion1'] > 0.1]
    markers[c] = cluster_df.index.tolist()[:3]

図でLeidenクラスタのデンドログラムを表示したいので事前に計算しておく。

In [58]:
sc.tl.dendrogram(adata, groupby='leiden_scVI', use_rep='X_scVI')

クラスタ間での遺伝子発現量の比較は、scanpy.pl.dotplotが見やすくて便利。上で抽出したマーカー遺伝子を指定して、表示する発現量の数値としてはadata.rawに格納したデータで描画する。

In [59]:
sc.pl.dotplot(
    adata,
    markers,
    groupby='leiden_scVI',
    dendrogram=True,
    color_map="Blues",
    swap_axes=True,
    use_raw=True,
    standard_scale='var')

クラスタ平均値ではなくクラスタごと細胞ごとの全部の（scVIノーマライズされた）発現量でヒートマップを描くこともできる。

In [60]:
sc.pl.heatmap(
    adata,
    markers,
    groupby='leiden_scVI',
    layer='scvi_normalized',
    standard_scale="var",
    dendrogram=True,
    figsize=(8, 12)
)

データの保存¶

In [61]:
adata.write(filename='./data/retinal.h5ad')
In [ ]: